home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 1.iso / HENSA / MATHS / PLPLOT / PLPLOT.ZIP / src / plot3d.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-07-20  |  31.9 KB  |  1,266 lines
  1. /* $Id: plot3d.c,v 1.16 1994/07/20 06:08:36 mjl Exp $
  2.  * $Log: plot3d.c,v $
  3.  * Revision 1.16  1994/07/20  06:08:36  mjl
  4.  * The new 3d functions moved elsewhere.
  5.  *
  6.  * Revision 1.15  1994/07/19  22:13:36  furnish
  7.  * Added new routine pl3poly() for drawing polygons in 3-space, with
  8.  * hidden surface removal.
  9.  *
  10.  * Revision 1.14  1994/07/18  20:33:30  mjl
  11.  * Some more cleaning up.
  12.  *
  13.  * Revision 1.13  1994/07/15  20:37:21  furnish
  14.  * Added routines pl3line and pl3poin for drawing lines and points in 3
  15.  * space.  Added a new example program, and dependency info to build it.
  16.  *
  17.  * Revision 1.12  1994/06/30  18:22:13  mjl
  18.  * All core source files: made another pass to eliminate warnings when using
  19.  * gcc -Wall.  Lots of cleaning up: got rid of includes of math.h or string.h
  20.  * (now included by plplot.h), and other minor changes.  Now each file has
  21.  * global access to the plstream pointer via extern; many accessor functions
  22.  * eliminated as a result.
  23.  *
  24.  * Revision 1.11  1994/05/13  22:56:35  mjl
  25.  * Fixed an old bug -- it was innocuous in the company of the allocation bugs
  26.  * fixed in the last update, which is why it wasn't discovered before now.
  27.  * Now plots correctly as well as frees all memory allocated.
  28.  *
  29.  * Revision 1.10  1994/04/08  12:34:21  mjl
  30.  * Fixed some cases of allocating memory and never freeing it.  Also changed
  31.  * some previously fatal errors into recoverable ones.
  32.  *
  33.  * Revision 1.9  1994/03/23  08:22:00  mjl
  34.  * Cruft elimination.
  35. */
  36.  
  37. /*    plot3d.c
  38.  
  39.     3d plot routines.
  40. */
  41.  
  42. #include "plplotP.h"
  43.  
  44. /* Internal constants */
  45.  
  46. #define  BINC      50        /* Block size for memory allocation */
  47.  
  48. static PLINT pl3mode = 0;    /* 0 3d solid; 1 mesh plot */
  49. static PLINT pl3upv = 1;    /* 1 update view; 0 no update */
  50.  
  51. static PLINT zbflg = 0, zbcol;
  52. static PLFLT zbtck;
  53.  
  54. static PLINT *oldhiview;
  55. static PLINT *oldloview;
  56. static PLINT *newhiview;
  57. static PLINT *newloview;
  58. static PLINT *u, *v;
  59.  
  60. static PLINT mhi, xxhi, newhisize;
  61. static PLINT mlo, xxlo, newlosize;
  62.  
  63. /* Prototypes for static functions */
  64. /* INDENT OFF */
  65.  
  66. static void plgrid3    (PLFLT);
  67. static void plnxtv    (PLINT *, PLINT *, PLINT, PLINT);
  68. static void plside3    (PLFLT *, PLFLT *, PLFLT **, PLINT, PLINT, PLINT);
  69. static void plt3zz    (PLINT, PLINT, PLINT, PLINT, 
  70.              PLINT, PLINT *, PLFLT *, PLFLT *, PLFLT **, 
  71.              PLINT, PLINT, PLINT *, PLINT *);
  72. static void plnxtvhi    (PLINT *, PLINT *, PLINT, PLINT);
  73. static void plnxtvlo    (PLINT *, PLINT *, PLINT, PLINT);
  74. static void savehipoint    (PLINT, PLINT);
  75. static void savelopoint    (PLINT, PLINT);
  76. static void swaphiview    (void);
  77. static void swaploview    (void);
  78. static void myexit    (char *);
  79. static void myabort    (char *);
  80. static void freework    (void);
  81. static int  plabv    (PLINT, PLINT, PLINT, PLINT, PLINT, PLINT);
  82. static void pl3cut    (PLINT, PLINT, PLINT, PLINT, PLINT, 
  83.                 PLINT, PLINT, PLINT, PLINT *, PLINT *);
  84. /* INDENT ON */
  85.  
  86. /*----------------------------------------------------------------------*\
  87.  * void plmesh(x, y, z, nx, ny, opt)
  88.  *
  89.  * Plots a mesh representation of the function z[x][y]. The x values
  90.  * are stored as x[0..nx-1], the y values as y[0..ny-1], and the
  91.  * z values are in the 2-d array z[][]. The integer "opt" specifies:
  92.  *
  93.  *  opt = 1:  Draw lines parallel to x-axis
  94.  *  opt = 2:  Draw lines parallel to y-axis
  95.  *  opt = 3:  Draw lines parallel to both axes
  96. \*----------------------------------------------------------------------*/
  97.  
  98. void
  99. c_plmesh(PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny, PLINT opt)
  100. {
  101.     pl3mode = 1;
  102.     plot3d(x, y, z, nx, ny, opt, 0);
  103.  
  104.     pl3mode = 0;
  105. }
  106.  
  107. /*----------------------------------------------------------------------*\
  108.  * void plot3d(x, y, z, nx, ny, opt, side)
  109.  *
  110.  * Plots a 3-d representation of the function z[x][y]. The x values
  111.  * are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
  112.  * values are in the 2-d array z[][]. The integer "opt" specifies:
  113.  *
  114.  *  opt = 1:  Draw lines parallel to x-axis
  115.  *  opt = 2:  Draw lines parallel to y-axis
  116.  *  opt = 3:  Draw lines parallel to both axes
  117. \*----------------------------------------------------------------------*/
  118.  
  119. void
  120. c_plot3d(PLFLT *x, PLFLT *y, PLFLT **z,
  121.      PLINT nx, PLINT ny, PLINT opt, PLINT side)
  122. {
  123.     PLINT color, font;
  124.     PLFLT cxx, cxy, cyx, cyy, cyz;
  125.     PLINT init, i, ix, iy;
  126.  
  127.     if (plsc->level < 3) {
  128.     myabort("plot3d: Please set up window first");
  129.     return;
  130.     }
  131.  
  132.     if (opt < 1 || opt > 3) {
  133.     myabort("plot3d: Bad option");
  134.     return;
  135.     }
  136.  
  137.     if (nx <= 0 || ny <= 0) {
  138.     myabort("plot3d: Bad array dimensions.");
  139.     return;
  140.     }
  141.  
  142. /* Check that points in x and in y are strictly increasing */
  143.  
  144.     for (i = 0; i < nx - 1; i++) {
  145.     if (x[i] >= x[i + 1]) {
  146.         myabort("plot3d: X array must be strictly increasing");
  147.         return;
  148.     }
  149.     }
  150.     for (i = 0; i < ny - 1; i++) {
  151.     if (y[i] >= y[i + 1]) {
  152.         myabort("plot3d: Y array must be strictly increasing");
  153.         return;
  154.     }
  155.     }
  156.  
  157. /* Allocate work arrays */
  158.  
  159.     u = (PLINT *) malloc((size_t) (2 * MAX(nx, ny) * sizeof(PLINT)));
  160.     v = (PLINT *) malloc((size_t) (2 * MAX(nx, ny) * sizeof(PLINT)));
  161.     if ( ! u || ! v)
  162.     myexit("plot3d: Out of memory.");
  163.  
  164.     plP_gw3wc(&cxx, &cxy, &cyx, &cyy, &cyz);
  165.     init = 1;
  166.  
  167. /* Call 3d line plotter.  Each viewing quadrant (perpendicular to x-y */
  168. /* plane) must be handled separately. */ 
  169.  
  170.     if (cxx >= 0.0 && cxy <= 0.0) {
  171.     if (opt == 2) 
  172.         plt3zz(1, ny, 1, -1, -opt, &init, x, y, z, nx, ny, u, v);
  173.     else {
  174.         for (iy = 2; iy <= ny; iy++)
  175.         plt3zz(1, iy, 1, -1, -opt, &init, x, y, z, nx, ny, u, v);
  176.     }
  177.     if (opt == 1)
  178.         plt3zz(1, ny, 1, -1, opt, &init, x, y, z, nx, ny, u, v);
  179.     else {
  180.         for (ix = 1; ix <= nx - 1; ix++)
  181.         plt3zz(ix, ny, 1, -1, opt, &init, x, y, z, nx, ny, u, v);
  182.     }
  183.     }
  184.  
  185.     else if (cxx <= 0.0 && cxy <= 0.0) {
  186.     if (opt == 1) 
  187.         plt3zz(nx, ny, -1, -1, opt, &init, x, y, z, nx, ny, u, v);
  188.     else {
  189.         for (ix = 2; ix <= nx; ix++) 
  190.         plt3zz(ix, ny, -1, -1, opt, &init, x, y, z, nx, ny, u, v);
  191.     }
  192.     if (opt == 2)
  193.         plt3zz(nx, ny, -1, -1, -opt, &init, x, y, z, nx, ny, u, v);
  194.     else {
  195.         for (iy = ny; iy >= 2; iy--)
  196.         plt3zz(nx, iy, -1, -1, -opt, &init, x, y, z, nx, ny, u, v);
  197.     }
  198.     }
  199.  
  200.     else if (cxx <= 0.0 && cxy >= 0.0) {
  201.     if (opt == 2) 
  202.         plt3zz(nx, 1, -1, 1, -opt, &init, x, y, z, nx, ny, u, v);
  203.     else {
  204.         for (iy = ny - 1; iy >= 1; iy--) 
  205.         plt3zz(nx, iy, -1, 1, -opt, &init, x, y, z, nx, ny, u, v);
  206.     }
  207.     if (opt == 1)
  208.         plt3zz(nx, 1, -1, 1, opt, &init, x, y, z, nx, ny, u, v);
  209.     else {
  210.         for (ix = nx; ix >= 2; ix--)
  211.         plt3zz(ix, 1, -1, 1, opt, &init, x, y, z, nx, ny, u, v);
  212.     }
  213.     }
  214.  
  215.     else if (cxx >= 0.0 && cxy >= 0.0) {
  216.     if (opt == 1) 
  217.         plt3zz(1, 1, 1, 1, opt, &init, x, y, z, nx, ny, u, v);
  218.     else {
  219.         for (ix = nx - 1; ix >= 1; ix--) 
  220.         plt3zz(ix, 1, 1, 1, opt, &init, x, y, z, nx, ny, u, v);
  221.     }
  222.     if (opt == 2)
  223.         plt3zz(1, 1, 1, 1, -opt, &init, x, y, z, nx, ny, u, v);
  224.     else {
  225.         for (iy = 1; iy <= ny - 1; iy++)
  226.         plt3zz(1, iy, 1, 1, -opt, &init, x, y, z, nx, ny, u, v);
  227.     }
  228.     }
  229.  
  230. /* Finish up by drawing sides, background grid (both are optional) */
  231.  
  232.     if (side)
  233.     plside3(x, y, z, nx, ny, opt);
  234.  
  235.     if (zbflg) {
  236.     plP_gatt(&font, &color);
  237.     plcol(zbcol);
  238.     plgrid3(zbtck);
  239.     plcol(color);
  240.     }
  241.  
  242.     freework();
  243. }
  244.  
  245. /*----------------------------------------------------------------------*\
  246.  * void plP_gzback()
  247.  *
  248.  * Get background parameters for 3d plot.
  249. \*----------------------------------------------------------------------*/
  250.  
  251. void
  252. plP_gzback(PLINT **zbf, PLINT **zbc, PLFLT **zbt)
  253. {
  254.     *zbf = &zbflg;
  255.     *zbc = &zbcol;
  256.     *zbt = &zbtck;
  257. }
  258.  
  259. /*----------------------------------------------------------------------*\
  260.  * void plt3zz()
  261.  *
  262.  * Draws the next zig-zag line for a 3-d plot. The data is stored in array
  263.  * z[][] as a function of x[] and y[]. This is plotted out starting at
  264.  * index (xstar0,ystar0).
  265.  *
  266.  * Depending on the state of "flg0", the sequence of data points sent to
  267.  * plnxtv is altered so as to allow cross-hatch plotting, or plotting
  268.  * parallel to either the x-axis or the y-axis.
  269. \*----------------------------------------------------------------------*/
  270.  
  271. static void
  272. plt3zz(PLINT xstar0, PLINT ystar0, PLINT dx, PLINT dy, PLINT flg0, PLINT *init,
  273.        PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny, PLINT *u, PLINT *v)
  274. {
  275.     PLINT flag;
  276.     PLINT n;
  277.     PLINT xstart, ystart;
  278.  
  279.     n = 0;
  280.     xstart = xstar0;
  281.     ystart = ystar0;
  282.     flag = flg0;
  283.  
  284.     while (1 <= xstart && xstart <= nx && 1 <= ystart && ystart <= ny) {
  285.     u[n] = plP_wcpcx(plP_w3wcx(
  286.             x[xstart - 1], y[ystart - 1], z[xstart - 1][ystart - 1]));
  287.  
  288.     v[n] = plP_wcpcy(plP_w3wcy(
  289.             x[xstart - 1], y[ystart - 1], z[xstart - 1][ystart - 1]));
  290.  
  291.     if (flag == -3) {
  292.         ystart = ystart + dy;
  293.         flag = -flag;
  294.     }
  295.     else if (flag == -2)
  296.         ystart = ystart + dy;
  297.  
  298.     else if (flag == -1) {
  299.         ystart = ystart + dy;
  300.         flag = -flag;
  301.     }
  302.     else if (flag == 1)
  303.         xstart = xstart + dx;
  304.  
  305.     else if (flag == 2) {
  306.         xstart = xstart + dx;
  307.         flag = -flag;
  308.     }
  309.     else if (flag == 3) {
  310.         xstart = xstart + dx;
  311.         flag = -flag;
  312.     }
  313.     n = n + 1;
  314.     }
  315.  
  316.     if (flag == 1 || flag == -2) {
  317.     if (flag == 1) {
  318.         xstart = xstart - dx;
  319.         ystart = ystart + dy;
  320.     }
  321.     else if (flag == -2) {
  322.         ystart = ystart - dy;
  323.         xstart = xstart + dx;
  324.     }
  325.  
  326.     if (1 <= xstart && xstart <= nx && 1 <= ystart && ystart <= ny) {
  327.         u[n] = plP_wcpcx(plP_w3wcx(
  328.             x[xstart - 1], y[ystart - 1], z[xstart - 1][ystart - 1]));
  329.  
  330.         v[n] = plP_wcpcy(plP_w3wcy(
  331.             x[xstart - 1], y[ystart - 1], z[xstart - 1][ystart - 1]));
  332.  
  333.         n = n + 1;
  334.     }
  335.     }
  336.     plnxtv(u, v, n, *init);
  337.     *init = 0;
  338. }
  339.  
  340. /*----------------------------------------------------------------------*\
  341.  * void plside3()
  342.  *
  343.  * This routine draws sides around the front of the 3d plot so that
  344.  * it does not appear to float.
  345. \*----------------------------------------------------------------------*/
  346.  
  347. static void
  348. plside3(PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny, PLINT opt)
  349. {
  350.     PLINT i;
  351.     PLFLT cxx, cxy, cyx, cyy, cyz;
  352.     PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
  353.     PLFLT tx, ty, ux, uy;
  354.  
  355.     plP_gw3wc(&cxx, &cxy, &cyx, &cyy, &cyz);
  356.     plP_gdom(&xmin, &xmax, &ymin, &ymax);
  357.     plP_grange(&zscale, &zmin, &zmax);
  358.  
  359.     if (cxx >= 0.0 && cxy <= 0.0) {
  360.     /* Get x, y coordinates of legs and plot */
  361.     if (opt != 1) {
  362.         for (i = 0; i < nx; i++) {
  363.         tx = plP_w3wcx(x[i], y[0], zmin);
  364.         ty = plP_w3wcy(x[i], y[0], zmin);
  365.         ux = plP_w3wcx(x[i], y[0], z[i][0]);
  366.         uy = plP_w3wcy(x[i], y[0], z[i][0]);
  367.         pljoin(tx, ty, ux, uy);
  368.         }
  369.     }
  370.  
  371.     if (opt != 2) {
  372.         for (i = 0; i < ny; i++) {
  373.         tx = plP_w3wcx(x[0], y[i], zmin);
  374.         ty = plP_w3wcy(x[0], y[i], zmin);
  375.         ux = plP_w3wcx(x[0], y[i], z[0][i]);
  376.         uy = plP_w3wcy(x[0], y[i], z[0][i]);
  377.         pljoin(tx, ty, ux, uy);
  378.         }
  379.     }
  380.     }
  381.     else if (cxx <= 0.0 && cxy <= 0.0) {
  382.     if (opt != 1) {
  383.         for (i = 0; i < nx; i++) {
  384.         tx = plP_w3wcx(x[i], y[ny - 1], zmin);
  385.         ty = plP_w3wcy(x[i], y[ny - 1], zmin);
  386.         ux = plP_w3wcx(x[i], y[ny - 1], z[i][ny - 1]);
  387.         uy = plP_w3wcy(x[i], y[ny - 1], z[i][ny - 1]);
  388.         pljoin(tx, ty, ux, uy);
  389.         }
  390.     }
  391.  
  392.     if (opt != 2) {
  393.         for (i = 0; i < ny; i++) {
  394.         tx = plP_w3wcx(x[0], y[i], zmin);
  395.         ty = plP_w3wcy(x[0], y[i], zmin);
  396.         ux = plP_w3wcx(x[0], y[i], z[0][i]);
  397.         uy = plP_w3wcy(x[0], y[i], z[0][i]);
  398.         pljoin(tx, ty, ux, uy);
  399.         }
  400.     }
  401.     }
  402.     else if (cxx <= 0.0 && cxy >= 0.0) {
  403.     if (opt != 1) {
  404.         for (i = 0; i < nx; i++) {
  405.         tx = plP_w3wcx(x[i], y[ny - 1], zmin);
  406.         ty = plP_w3wcy(x[i], y[ny - 1], zmin);
  407.         ux = plP_w3wcx(x[i], y[ny - 1], z[i][ny - 1]);
  408.         uy = plP_w3wcy(x[i], y[ny - 1], z[i][ny - 1]);
  409.         pljoin(tx, ty, ux, uy);
  410.         }
  411.     }
  412.  
  413.     if (opt != 2) {
  414.         for (i = 0; i < ny; i++) {
  415.         tx = plP_w3wcx(x[nx - 1], y[i], zmin);
  416.         ty = plP_w3wcy(x[nx - 1], y[i], zmin);
  417.         ux = plP_w3wcx(x[nx - 1], y[i], z[nx - 1][i]);
  418.         uy = plP_w3wcy(x[nx - 1], y[i], z[nx - 1][i]);
  419.         pljoin(tx, ty, ux, uy);
  420.         }
  421.     }
  422.     }
  423.     else if (cxx >= 0.0 && cxy >= 0.0) {
  424.     if (opt != 1) {
  425.         for (i = 0; i < nx; i++) {
  426.         tx = plP_w3wcx(x[i], y[0], zmin);
  427.         ty = plP_w3wcy(x[i], y[0], zmin);
  428.         ux = plP_w3wcx(x[i], y[0], z[i][0]);
  429.         uy = plP_w3wcy(x[i], y[0], z[i][0]);
  430.         pljoin(tx, ty, ux, uy);
  431.         }
  432.     }
  433.  
  434.     if (opt != 2) {
  435.         for (i = 0; i < ny; i++) {
  436.         tx = plP_w3wcx(x[nx - 1], y[i], zmin);
  437.         ty = plP_w3wcy(x[nx - 1], y[i], zmin);
  438.         ux = plP_w3wcx(x[nx - 1], y[i], z[nx - 1][i]);
  439.         uy = plP_w3wcy(x[nx - 1], y[i], z[nx - 1][i]);
  440.         pljoin(tx, ty, ux, uy);
  441.         }
  442.     }
  443.     }
  444. }
  445.  
  446. /*----------------------------------------------------------------------*\
  447.  * void plgrid3()
  448.  *
  449.  * This routine draws a grid around the back side of the 3d plot with
  450.  * hidden line removal.
  451. \*----------------------------------------------------------------------*/
  452.  
  453. static void
  454. plgrid3(PLFLT tick)
  455. {
  456.     PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
  457.     PLFLT cxx, cxy, cyx, cyy, cyz;
  458.     PLINT u[3], v[3];
  459.     PLINT nsub, prec, mode, digmax, digits, scale;
  460.     PLFLT tp;
  461.  
  462.     plP_gw3wc(&cxx, &cxy, &cyx, &cyy, &cyz);
  463.     plP_gdom(&xmin, &xmax, &ymin, &ymax);
  464.     plP_grange(&zscale, &zmin, &zmax);
  465.  
  466.     plgzax(&digmax, &digits);
  467.     nsub = 0;
  468.  
  469.     pldtik(zmin, zmax, &tick, &nsub, &mode, &prec, digmax, &scale);
  470.     tp = tick * floor(zmin / tick) + tick;
  471.     pl3upv = 0;
  472.  
  473.     if (cxx >= 0.0 && cxy <= 0.0) {
  474.     while (tp <= zmax) {
  475.         u[0] = plP_wcpcx(plP_w3wcx(xmin, ymax, tp));
  476.         v[0] = plP_wcpcy(plP_w3wcy(xmin, ymax, tp));
  477.         u[1] = plP_wcpcx(plP_w3wcx(xmax, ymax, tp));
  478.         v[1] = plP_wcpcy(plP_w3wcy(xmax, ymax, tp));
  479.         u[2] = plP_wcpcx(plP_w3wcx(xmax, ymin, tp));
  480.         v[2] = plP_wcpcy(plP_w3wcy(xmax, ymin, tp));
  481.         plnxtv(u, v, 3, 0);
  482.  
  483.         tp += tick;
  484.     }
  485.     u[0] = plP_wcpcx(plP_w3wcx(xmax, ymax, zmin));
  486.     v[0] = plP_wcpcy(plP_w3wcy(xmax, ymax, zmin));
  487.     u[1] = plP_wcpcx(plP_w3wcx(xmax, ymax, zmax));
  488.     v[1] = plP_wcpcy(plP_w3wcy(xmax, ymax, zmax));
  489.     plnxtv(u, v, 2, 0);
  490.     }
  491.     else if (cxx <= 0.0 && cxy <= 0.0) {
  492.     while (tp <= zmax) {
  493.         u[0] = plP_wcpcx(plP_w3wcx(xmax, ymax, tp));
  494.         v[0] = plP_wcpcy(plP_w3wcy(xmax, ymax, tp));
  495.         u[1] = plP_wcpcx(plP_w3wcx(xmax, ymin, tp));
  496.         v[1] = plP_wcpcy(plP_w3wcy(xmax, ymin, tp));
  497.         u[2] = plP_wcpcx(plP_w3wcx(xmin, ymin, tp));
  498.         v[2] = plP_wcpcy(plP_w3wcy(xmin, ymin, tp));
  499.         plnxtv(u, v, 3, 0);
  500.  
  501.         tp += tick;
  502.     }
  503.     u[0] = plP_wcpcx(plP_w3wcx(xmax, ymin, zmin));
  504.     v[0] = plP_wcpcy(plP_w3wcy(xmax, ymin, zmin));
  505.     u[1] = plP_wcpcx(plP_w3wcx(xmax, ymin, zmax));
  506.     v[1] = plP_wcpcy(plP_w3wcy(xmax, ymin, zmax));
  507.     plnxtv(u, v, 2, 0);
  508.     }
  509.     else if (cxx <= 0.0 && cxy >= 0.0) {
  510.     while (tp <= zmax) {
  511.         u[0] = plP_wcpcx(plP_w3wcx(xmax, ymin, tp));
  512.         v[0] = plP_wcpcy(plP_w3wcy(xmax, ymin, tp));
  513.         u[1] = plP_wcpcx(plP_w3wcx(xmin, ymin, tp));
  514.         v[1] = plP_wcpcy(plP_w3wcy(xmin, ymin, tp));
  515.         u[2] = plP_wcpcx(plP_w3wcx(xmin, ymax, tp));
  516.         v[2] = plP_wcpcy(plP_w3wcy(xmin, ymax, tp));
  517.         plnxtv(u, v, 3, 0);
  518.  
  519.         tp += tick;
  520.     }
  521.     u[0] = plP_wcpcx(plP_w3wcx(xmin, ymin, zmin));
  522.     v[0] = plP_wcpcy(plP_w3wcy(xmin, ymin, zmin));
  523.     u[1] = plP_wcpcx(plP_w3wcx(xmin, ymin, zmax));
  524.     v[1] = plP_wcpcy(plP_w3wcy(xmin, ymin, zmax));
  525.     plnxtv(u, v, 2, 0);
  526.     }
  527.     else if (cxx >= 0.0 && cxy >= 0.0) {
  528.     while (tp <= zmax) {
  529.         u[0] = plP_wcpcx(plP_w3wcx(xmin, ymin, tp));
  530.         v[0] = plP_wcpcy(plP_w3wcy(xmin, ymin, tp));
  531.         u[1] = plP_wcpcx(plP_w3wcx(xmin, ymax, tp));
  532.         v[1] = plP_wcpcy(plP_w3wcy(xmin, ymax, tp));
  533.         u[2] = plP_wcpcx(plP_w3wcx(xmax, ymax, tp));
  534.         v[2] = plP_wcpcy(plP_w3wcy(xmax, ymax, tp));
  535.         plnxtv(u, v, 3, 0);
  536.  
  537.         tp += tick;
  538.     }
  539.     u[0] = plP_wcpcx(plP_w3wcx(xmin, ymax, zmin));
  540.     v[0] = plP_wcpcy(plP_w3wcy(xmin, ymax, zmin));
  541.     u[1] = plP_wcpcx(plP_w3wcx(xmin, ymax, zmax));
  542.     v[1] = plP_wcpcy(plP_w3wcy(xmin, ymax, zmax));
  543.     plnxtv(u, v, 2, 0);
  544.     }
  545.     pl3upv = 1;
  546. }
  547.  
  548. /*----------------------------------------------------------------------*\
  549.  * void plnxtv()
  550.  *
  551.  * Draw the next view of a 3-d plot. The physical coordinates of the
  552.  * points for the next view are placed in the n points of arrays u and
  553.  * v. The silhouette found so far is stored in the heap as a set of m peak
  554.  * points.
  555.  *
  556.  * These routines dynamically allocate memory for hidden line removal.
  557.  * Memory is allocated in blocks of 2*BINC*sizeof(PLINT) bytes.  Large
  558.  * values of BINC give better performance but also allocate more memory
  559.  * than is needed. If your 3D plots are very "spiky" or you are working
  560.  * with very large matrices then you will probably want to increase BINC.
  561. \*----------------------------------------------------------------------*/
  562.  
  563. static void
  564. plnxtv(PLINT *u, PLINT *v, PLINT n, PLINT init)
  565. {
  566.     plnxtvhi(u, v, n, init);
  567.  
  568.     if (pl3mode)
  569.     plnxtvlo(u, v, n, init);
  570. }
  571.  
  572. /*----------------------------------------------------------------------*\
  573.  * void plnxtvhi()
  574.  *
  575.  * Draw the top side of the 3-d plot.
  576. \*----------------------------------------------------------------------*/
  577.  
  578. static void
  579. plnxtvhi(PLINT *u, PLINT *v, PLINT n, PLINT init)
  580. {
  581.     PLINT i, j, first;
  582.     PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
  583.     PLINT su1, su2, sv1, sv2;
  584.     PLINT cx, cy, px, py;
  585.     PLINT seg, ptold, lstold = 0, pthi, pnewhi, newhi, change, ochange = 0;
  586.  
  587.     first = 1;
  588.     pnewhi = 0;
  589.  
  590. /*
  591.  * For the initial set of points, just display them and store them as the
  592.  * peak points.
  593.  */
  594.     if (init == 1) {
  595.  
  596.     oldhiview = (PLINT *) malloc((size_t) (2 * n * sizeof(PLINT)));
  597.     if ( ! oldhiview)
  598.         myexit("plnxtvhi: Out of memory.");
  599.  
  600.     plP_movphy(u[0], v[0]);
  601.     oldhiview[0] = u[0];
  602.     oldhiview[1] = v[0];
  603.     for (i = 1; i < n; i++) {
  604.         plP_draphy(u[i], v[i]);
  605.         oldhiview[2 * i] = u[i];
  606.         oldhiview[2 * i + 1] = v[i];
  607.     }
  608.     mhi = n;
  609.     return;
  610.     }
  611.  
  612. /*
  613.  * Otherwise, we need to consider hidden-line removal problem. We scan
  614.  * over the points in both the old (i.e. oldhiview[]) and new (i.e. u[]
  615.  * and v[]) arrays in order of increasing x coordinate.  At each stage, we
  616.  * find the line segment in the other array (if one exists) that straddles
  617.  * the x coordinate of the point. We have to determine if the point lies
  618.  * above or below the line segment, and to check if the below/above status
  619.  * has changed since the last point.
  620.  *
  621.  * If pl3upv = 0 we do not update the view, this is useful for drawing
  622.  * lines on the graph after we are done plotting points.  Hidden line
  623.  * removal is still done, but the view is not updated.
  624.  */
  625.     xxhi = 0;
  626.     i = 0;
  627.     j = 0;
  628.     if (pl3upv != 0) {
  629.     newhisize = 2 * (mhi + BINC);
  630.     if (newhiview != NULL) {
  631.         newhiview = 
  632.         (PLINT *) realloc((void *) newhiview,
  633.                   (size_t) (newhisize * sizeof(PLINT)));
  634.     }
  635.     else {
  636.         newhiview = 
  637.         (PLINT *) malloc((size_t) (newhisize * sizeof(PLINT)));
  638.     }
  639.     if ( ! newhiview)
  640.         myexit("plnxtvhi: Out of memory.");
  641.     }
  642.  
  643. /*
  644.  * (oldhiview[2*i], oldhiview[2*i]) is the i'th point in the old array
  645.  * (u[j], v[j]) is the j'th point in the new array
  646.  */
  647.     while (i < mhi || j < n) {
  648.  
  649. /*
  650.  * The coordinates of the point under consideration are (px,py).  The line
  651.  * segment joins (sx1,sy1) to (sx2,sy2).  "ptold" is true if the point
  652.  * lies in the old array. We set it by comparing the x coordinates of the
  653.  * i'th old point and the j'th new point, being careful if we have fallen
  654.  * past the edges. Having found the point, load up the point and segment
  655.  * coordinates appropriately.
  656.  */
  657.     ptold = ((oldhiview[2 * i] < u[j] && i < mhi) || j >= n);
  658.     if (ptold) {
  659.         px = oldhiview[2 * i];
  660.         py = oldhiview[2 * i + 1];
  661.         seg = j > 0 && j < n;
  662.         if (seg) {
  663.         sx1 = u[j - 1];
  664.         sy1 = v[j - 1];
  665.         sx2 = u[j];
  666.         sy2 = v[j];
  667.         }
  668.     }
  669.     else {
  670.         px = u[j];
  671.         py = v[j];
  672.         seg = i > 0 && i < mhi;
  673.         if (seg) {
  674.         sx1 = oldhiview[2 * (i - 1)];
  675.         sy1 = oldhiview[2 * (i - 1) + 1];
  676.         sx2 = oldhiview[2 * i];
  677.         sy2 = oldhiview[2 * i + 1];
  678.         }
  679.     }
  680.  
  681. /*
  682.  * Now determine if the point is higher than the segment, using the
  683.  * logical function "above". We also need to know if it is the old view or
  684.  * the new view that is higher. "newhi" is set true if the new view is
  685.  * higher than the old.
  686.  */
  687.     if (seg)
  688.         pthi = plabv(px, py, sx1, sy1, sx2, sy2);
  689.     else
  690.         pthi = 1;
  691.  
  692.     newhi = (ptold && !pthi) || (!ptold && pthi);
  693.     change = (newhi && !pnewhi) || (!newhi && pnewhi);
  694.  
  695. /*
  696.  * There is a new intersection point to put in the peak array if the state
  697.  * of "newhi" changes.
  698.  */
  699.     if (first) {
  700.         plP_movphy(px, py);
  701.         first = 0;
  702.         lstold = ptold;
  703.         savehipoint(px, py);
  704.         pthi = 0;
  705.         ochange = 0;
  706.     }
  707.     else if (change) {
  708.  
  709. /*
  710.  * Take care of special cases at end of arrays.  If pl3upv is 0 the
  711.  * endpoints are not connected to the old view.
  712.  */
  713.         if (pl3upv == 0 && ((!ptold && j == 0) || (ptold && i == 0))) {
  714.         plP_movphy(px, py);
  715.         lstold = ptold;
  716.         pthi = 0;
  717.         ochange = 0;
  718.         }
  719.         else if (pl3upv == 0 &&
  720.              (( ! ptold && i >= mhi) || (ptold && j >= n))) {
  721.         plP_movphy(px, py);
  722.         lstold = ptold;
  723.         pthi = 0;
  724.         ochange = 0;
  725.         }
  726.  
  727. /*
  728.  * If pl3upv is not 0 then we do want to connect the current line with the
  729.  * previous view at the endpoints.  Also find intersection point with old
  730.  * view.
  731.  */
  732.         else {
  733.         if (i == 0) {
  734.             sx1 = oldhiview[0];
  735.             sy1 = -1;
  736.             sx2 = oldhiview[0];
  737.             sy2 = oldhiview[1];
  738.         }
  739.         else if (i >= mhi) {
  740.             sx1 = oldhiview[2 * (mhi - 1)];
  741.             sy1 = oldhiview[2 * (mhi - 1) + 1];
  742.             sx2 = oldhiview[2 * (mhi - 1)];
  743.             sy2 = -1;
  744.         }
  745.         else {
  746.             sx1 = oldhiview[2 * (i - 1)];
  747.             sy1 = oldhiview[2 * (i - 1) + 1];
  748.             sx2 = oldhiview[2 * i];
  749.             sy2 = oldhiview[2 * i + 1];
  750.         }
  751.  
  752.         if (j == 0) {
  753.             su1 = u[0];
  754.             sv1 = -1;
  755.             su2 = u[0];
  756.             sv2 = v[0];
  757.         }
  758.         else if (j >= n) {
  759.             su1 = u[n - 1];
  760.             sv1 = v[n - 1];
  761.             su2 = u[n];
  762.             sv2 = -1;
  763.         }
  764.         else {
  765.             su1 = u[j - 1];
  766.             sv1 = v[j - 1];
  767.             su2 = u[j];
  768.             sv2 = v[j];
  769.         }
  770.  
  771. /* Determine the intersection */
  772.  
  773.         pl3cut(sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy);
  774.         if (cx == px && cy == py) {
  775.             if (lstold && !ochange)
  776.             plP_movphy(px, py);
  777.             else
  778.             plP_draphy(px, py);
  779.  
  780.             savehipoint(px, py);
  781.             lstold = 1;
  782.             pthi = 0;
  783.         }
  784.         else {
  785.             if (lstold && !ochange)
  786.             plP_movphy(cx, cy);
  787.             else
  788.             plP_draphy(cx, cy);
  789.  
  790.             lstold = 1;
  791.             savehipoint(cx, cy);
  792.         }
  793.         ochange = 1;
  794.         }
  795.     }
  796.  
  797. /* If point is high then draw plot to point and update view. */
  798.  
  799.     if (pthi) {
  800.         if (lstold && ptold)
  801.         plP_movphy(px, py);
  802.         else
  803.         plP_draphy(px, py);
  804.  
  805.         savehipoint(px, py);
  806.         lstold = ptold;
  807.         ochange = 0;
  808.     }
  809.  
  810.     pnewhi = newhi;
  811.  
  812.     if (ptold)
  813.         i = i + 1;
  814.     else
  815.         j = j + 1;
  816.     }
  817.  
  818. /* Set oldhiview */
  819.  
  820.     swaphiview();
  821. }
  822.  
  823. /*----------------------------------------------------------------------*\
  824.  * void plnxtvlo()
  825.  *
  826.  * Draw the bottom side of the 3-d plot.
  827. \*----------------------------------------------------------------------*/
  828.  
  829. static void
  830. plnxtvlo(PLINT *u, PLINT *v, PLINT n, PLINT init)
  831. {
  832.     PLINT i, j, first;
  833.     PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
  834.     PLINT su1, su2, sv1, sv2;
  835.     PLINT cx, cy, px, py;
  836.     PLINT seg, ptold, lstold = 0, ptlo, pnewlo, newlo, change, ochange = 0;
  837.  
  838.     first = 1;
  839.     pnewlo = 0;
  840.  
  841. /*
  842.  * For the initial set of points, just display them and store them as the
  843.  * peak points.
  844.  */
  845.     if (init == 1) {
  846.  
  847.     oldloview = (PLINT *) malloc((size_t) (2 * n * sizeof(PLINT)));
  848.     if ( ! oldloview)
  849.         myexit("plnxtvlo: Out of memory.");
  850.  
  851.     plP_movphy(u[0], v[0]);
  852.     oldloview[0] = u[0];
  853.     oldloview[1] = v[0];
  854.     for (i = 1; i < n; i++) {
  855.         plP_draphy(u[i], v[i]);
  856.         oldloview[2 * i] = u[i];
  857.         oldloview[2 * i + 1] = v[i];
  858.     }
  859.     mlo = n;
  860.     return;
  861.     }
  862.  
  863. /*
  864.  * Otherwise, we need to consider hidden-line removal problem. We scan
  865.  * over the points in both the old (i.e. oldloview[]) and new (i.e. u[]
  866.  * and v[]) arrays in order of increasing x coordinate.  At each stage, we
  867.  * find the line segment in the other array (if one exists) that straddles
  868.  * the x coordinate of the point. We have to determine if the point lies
  869.  * above or below the line segment, and to check if the below/above status
  870.  * has changed since the last point.
  871.  *
  872.  * If pl3upv = 0 we do not update the view, this is useful for drawing
  873.  * lines on the graph after we are done plotting points.  Hidden line
  874.  * removal is still done, but the view is not updated.
  875.  */
  876.     xxlo = 0;
  877.     i = 0;
  878.     j = 0;
  879.     if (pl3upv != 0) {
  880.     newlosize = 2 * (mlo + BINC);
  881.     if (newloview != NULL) {
  882.         newloview = 
  883.         (PLINT *) realloc((void *) newloview,
  884.                   (size_t) (newlosize * sizeof(PLINT)));
  885.     }
  886.     else {
  887.         newloview = 
  888.         (PLINT *) malloc((size_t) (newlosize * sizeof(PLINT)));
  889.     }
  890.     if ( ! newloview)
  891.         myexit("plnxtvlo: Out of memory.");
  892.     }
  893.  
  894. /*
  895.  * (oldloview[2*i], oldloview[2*i]) is the i'th point in the old array
  896.  * (u[j], v[j]) is the j'th point in the new array.
  897.  */
  898.     while (i < mlo || j < n) {
  899.  
  900. /*
  901.  * The coordinates of the point under consideration are (px,py).  The line
  902.  * segment joins (sx1,sy1) to (sx2,sy2).  "ptold" is true if the point
  903.  * lies in the old array. We set it by comparing the x coordinates of the
  904.  * i'th old point and the j'th new point, being careful if we have fallen
  905.  * past the edges. Having found the point, load up the point and segment
  906.  * coordinates appropriately.
  907.  */
  908.     ptold = ((oldloview[2 * i] < u[j] && i < mlo) || j >= n);
  909.     if (ptold) {
  910.         px = oldloview[2 * i];
  911.         py = oldloview[2 * i + 1];
  912.         seg = j > 0 && j < n;
  913.         if (seg) {
  914.         sx1 = u[j - 1];
  915.         sy1 = v[j - 1];
  916.         sx2 = u[j];
  917.         sy2 = v[j];
  918.         }
  919.     }
  920.     else {
  921.         px = u[j];
  922.         py = v[j];
  923.         seg = i > 0 && i < mlo;
  924.         if (seg) {
  925.         sx1 = oldloview[2 * (i - 1)];
  926.         sy1 = oldloview[2 * (i - 1) + 1];
  927.         sx2 = oldloview[2 * i];
  928.         sy2 = oldloview[2 * i + 1];
  929.         }
  930.     }
  931.  
  932. /*
  933.  * Now determine if the point is lower than the segment, using the logical
  934.  * function "above". We also need to know if it is the old view or the new
  935.  * view that is lower. "newlo" is set true if the new view is lower than
  936.  * the old.
  937.  */
  938.     if (seg)
  939.         ptlo = !plabv(px, py, sx1, sy1, sx2, sy2);
  940.     else
  941.         ptlo = 1;
  942.  
  943.     newlo = (ptold && !ptlo) || (!ptold && ptlo);
  944.     change = (newlo && !pnewlo) || (!newlo && pnewlo);
  945.  
  946. /*
  947.  * There is a new intersection point to put in the peak array if the state
  948.  * of "newlo" changes.
  949.  */
  950.     if (first) {
  951.         plP_movphy(px, py);
  952.         first = 0;
  953.         lstold = ptold;
  954.         savelopoint(px, py);
  955.         ptlo = 0;
  956.         ochange = 0;
  957.     }
  958.     else if (change) {
  959.  
  960. /*
  961.  * Take care of special cases at end of arrays.  If pl3upv is 0 the
  962.  * endpoints are not connected to the old view.
  963.  */
  964.         if (pl3upv == 0 && ((!ptold && j == 0) || (ptold && i == 0))) {
  965.         plP_movphy(px, py);
  966.         lstold = ptold;
  967.         ptlo = 0;
  968.         ochange = 0;
  969.         }
  970.         else if (pl3upv == 0 &&
  971.              (( ! ptold && i >= mlo) || (ptold && j >= n))) {
  972.         plP_movphy(px, py);
  973.         lstold = ptold;
  974.         ptlo = 0;
  975.         ochange = 0;
  976.         }
  977.  
  978. /*
  979.  * If pl3upv is not 0 then we do want to connect the current line with the
  980.  * previous view at the endpoints.  Also find intersection point with old
  981.  * view.
  982.  */
  983.         else {
  984.         if (i == 0) {
  985.             sx1 = oldloview[0];
  986.             sy1 = 100000;
  987.             sx2 = oldloview[0];
  988.             sy2 = oldloview[1];
  989.         }
  990.         else if (i >= mlo) {
  991.             sx1 = oldloview[2 * (mlo - 1)];
  992.             sy1 = oldloview[2 * (mlo - 1) + 1];
  993.             sx2 = oldloview[2 * (mlo - 1)];
  994.             sy2 = 100000;
  995.         }
  996.         else {
  997.             sx1 = oldloview[2 * (i - 1)];
  998.             sy1 = oldloview[2 * (i - 1) + 1];
  999.             sx2 = oldloview[2 * i];
  1000.             sy2 = oldloview[2 * i + 1];
  1001.         }
  1002.  
  1003.         if (j == 0) {
  1004.             su1 = u[0];
  1005.             sv1 = 100000;
  1006.             su2 = u[0];
  1007.             sv2 = v[0];
  1008.         }
  1009.         else if (j >= n) {
  1010.             su1 = u[n - 1];
  1011.             sv1 = v[n - 1];
  1012.             su2 = u[n];
  1013.             sv2 = 100000;
  1014.         }
  1015.         else {
  1016.             su1 = u[j - 1];
  1017.             sv1 = v[j - 1];
  1018.             su2 = u[j];
  1019.             sv2 = v[j];
  1020.         }
  1021.  
  1022. /* Determine the intersection */
  1023.  
  1024.         pl3cut(sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy);
  1025.         if (cx == px && cy == py) {
  1026.             if (lstold && !ochange)
  1027.             plP_movphy(px, py);
  1028.             else
  1029.             plP_draphy(px, py);
  1030.  
  1031.             savelopoint(px, py);
  1032.             lstold = 1;
  1033.             ptlo = 0;
  1034.         }
  1035.         else {
  1036.             if (lstold && !ochange)
  1037.             plP_movphy(cx, cy);
  1038.             else
  1039.             plP_draphy(cx, cy);
  1040.  
  1041.             lstold = 1;
  1042.             savelopoint(cx, cy);
  1043.         }
  1044.         ochange = 1;
  1045.         }
  1046.     }
  1047.  
  1048. /* If point is low then draw plot to point and update view. */
  1049.  
  1050.     if (ptlo) {
  1051.         if (lstold && ptold)
  1052.         plP_movphy(px, py);
  1053.         else
  1054.         plP_draphy(px, py);
  1055.  
  1056.         savelopoint(px, py);
  1057.         lstold = ptold;
  1058.         ochange = 0;
  1059.     }
  1060.  
  1061.     pnewlo = newlo;
  1062.  
  1063.     if (ptold)
  1064.         i = i + 1;
  1065.     else
  1066.         j = j + 1;
  1067.     }
  1068.  
  1069. /* Set oldloview */
  1070.  
  1071.     swaploview();
  1072. }
  1073.  
  1074. /*----------------------------------------------------------------------*\
  1075.  * savehipoint
  1076.  * savelopoint
  1077.  *
  1078.  * Add a point to the list of currently visible peaks/valleys, when
  1079.  * drawing the top/bottom surface, respectively.
  1080. \*----------------------------------------------------------------------*/
  1081.  
  1082. static void
  1083. savehipoint(PLINT px, PLINT py)
  1084. {
  1085.     if (pl3upv == 0)
  1086.     return;
  1087.  
  1088.     if (xxhi >= newhisize) {    /* allocate additional space */
  1089.     newhisize += 2 * BINC;
  1090.     newhiview = (PLINT *) realloc((void *) newhiview,
  1091.                       (size_t) (newhisize * sizeof(PLINT)));
  1092.     if ( ! newhiview) 
  1093.         myexit("savehipoint: Out of memory.");
  1094.     }
  1095.  
  1096.     newhiview[xxhi] = px;
  1097.     xxhi++;
  1098.     newhiview[xxhi] = py;
  1099.     xxhi++;
  1100. }
  1101.  
  1102. static void
  1103. savelopoint(PLINT px, PLINT py)
  1104. {
  1105.     if (pl3upv == 0)
  1106.     return;
  1107.  
  1108.     if (xxlo >= newlosize) {    /* allocate additional space */
  1109.     newlosize += 2 * BINC;
  1110.     newloview = (PLINT *) realloc((void *) newloview,
  1111.                       (size_t) (newlosize * sizeof(PLINT)));
  1112.     if ( ! newloview)
  1113.         myexit("savelopoint: Out of memory.");
  1114.     }
  1115.  
  1116.     newloview[xxlo] = px;
  1117.     xxlo++;
  1118.     newloview[xxlo] = py;
  1119.     xxlo++;
  1120. }
  1121.  
  1122. /*----------------------------------------------------------------------*\
  1123.  * swaphiview
  1124.  * swaploview
  1125.  *
  1126.  * Swaps the top/bottom views.  Need to do a real swap so that the 
  1127.  * memory cleanup routine really frees everything (and only once).
  1128. \*----------------------------------------------------------------------*/
  1129.  
  1130. static void
  1131. swaphiview(void)
  1132. {
  1133.     PLINT *tmp;
  1134.  
  1135.     if (pl3upv != 0) {
  1136.     mhi = xxhi / 2;
  1137.     tmp = oldhiview;
  1138.     oldhiview = newhiview;
  1139.     newhiview = tmp;
  1140.     }
  1141. }
  1142.  
  1143. static void
  1144. swaploview(void)
  1145. {
  1146.     PLINT *tmp;
  1147.  
  1148.     if (pl3upv != 0) {
  1149.     mlo = xxlo / 2;
  1150.     tmp = oldloview;
  1151.     oldloview = newloview;
  1152.     newloview = tmp;
  1153.     }
  1154. }
  1155.  
  1156. /*----------------------------------------------------------------------*\
  1157.  * freework
  1158.  *
  1159.  * Frees memory associated with work arrays
  1160. \*----------------------------------------------------------------------*/
  1161.  
  1162. static void
  1163. freework(void)
  1164. {
  1165.     free_mem(oldhiview);
  1166.     free_mem(oldloview);
  1167.     free_mem(newhiview);
  1168.     free_mem(newloview);
  1169.     free_mem(v);
  1170.     free_mem(u);
  1171. }
  1172.  
  1173. /*----------------------------------------------------------------------*\
  1174.  * myexit
  1175.  *
  1176.  * Calls plexit, cleaning up first.
  1177. \*----------------------------------------------------------------------*/
  1178.  
  1179. static void
  1180. myexit(char *msg)
  1181. {
  1182.     freework();
  1183.     plexit(msg);
  1184. }
  1185.  
  1186. /*----------------------------------------------------------------------*\
  1187.  * myabort
  1188.  *
  1189.  * Calls plabort, cleaning up first.
  1190.  * Caller should return to the user program.
  1191. \*----------------------------------------------------------------------*/
  1192.  
  1193. static void
  1194. myabort(char *msg)
  1195. {
  1196.     freework();
  1197.     plabort(msg);
  1198. }
  1199.  
  1200. /*----------------------------------------------------------------------*\
  1201.  * int plabv()
  1202.  *
  1203.  * Determines if point (px,py) lies above the line joining (sx1,sy1) to
  1204.  * (sx2,sy2). It only works correctly if sx1 <= px <= sx2.
  1205. \*----------------------------------------------------------------------*/
  1206.  
  1207. static int
  1208. plabv(PLINT px, PLINT py, PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2)
  1209. {
  1210.     PLINT above;
  1211.  
  1212.     if (py >= sy1 && py >= sy2)
  1213.     above = 1;
  1214.     else if (py < sy1 && py < sy2)
  1215.     above = 0;
  1216.     else if ((PLFLT) (sx2 - sx1) * (py - sy1) >
  1217.          (PLFLT) (px - sx1) * (sy2 - sy1))
  1218.     above = 1;
  1219.     else
  1220.     above = 0;
  1221.  
  1222.     return ((PLINT) above);
  1223. }
  1224.  
  1225. /*----------------------------------------------------------------------*\
  1226.  * void pl3cut()
  1227.  *
  1228.  * Determines the point of intersection (cx,cy) between the line joining
  1229.  * (sx1,sy1) to (sx2,sy2) and the line joining (su1,sv1) to (su2,sv2).
  1230. \*----------------------------------------------------------------------*/
  1231.  
  1232. static void
  1233. pl3cut(PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2,
  1234.        PLINT su1, PLINT sv1, PLINT su2, PLINT sv2, PLINT *cx, PLINT *cy)
  1235. {
  1236.     PLINT x21, y21, u21, v21, yv1, xu1, a, b;
  1237.     PLFLT fa, fb;
  1238.  
  1239.     x21 = sx2 - sx1;
  1240.     y21 = sy2 - sy1;
  1241.     u21 = su2 - su1;
  1242.     v21 = sv2 - sv1;
  1243.     yv1 = sy1 - sv1;
  1244.     xu1 = sx1 - su1;
  1245.  
  1246.     a = x21 * v21 - y21 * u21;
  1247.     fa = (PLFLT) a;
  1248.     if (a == 0) {
  1249.     if (sx2 < su2) {
  1250.         *cx = sx2;
  1251.         *cy = sy2;
  1252.     }
  1253.     else {
  1254.         *cx = su2;
  1255.         *cy = sv2;
  1256.     }
  1257.     return;
  1258.     }
  1259.     else {
  1260.     b = yv1 * u21 - xu1 * v21;
  1261.     fb = (PLFLT) b;
  1262.     *cx = (PLINT) (sx1 + (fb * x21) / fa + .5);
  1263.     *cy = (PLINT) (sy1 + (fb * y21) / fa + .5);
  1264.     }
  1265. }
  1266.